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The universal equations describing collective oscillations of the multidomain patterns of small 
period in an arbitrary d- dimensional reaction-diffusion system of the activator-inhibitor type are 
asymptotically derived. It is shown that not far from the instability leading to the formation 
of the pulsating multidomain pattern the oscillations of different domains synchronize. In one 
dimension standing and traveling waves of the oscillation phase are realized. In addition to these, 
in two dimensions target and spiral waves of the oscillation phase, as well as spatio-temporal chaos 
of domain oscillations, are feasible. Further inside the unstable region the collective oscillations 
break down, so the pulsating multidomain pattern transforms into an irregular pulsating pattern, 
the uniform self-oscillations, or turbulence. The parameter regions where these effects occur are 
analyzed. The effects of the pattern's disorder are also studied. The conclusions of the analysis 
are supported by the numerical simulations of a concrete model. The obtained results explain the 
dynamics of Turing patterns observed in the experiments on chlorite-iodide-malonic acid reaction. 
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I. INTRODUCTION 

Complex dynamic patterns, such as traveling waves, 
breathers, or spatio-temporal chaos, are encountered in 
the variety of noncquilibrium physical, chemical, and bi- 
ological systems Jl|--pT|. These systems include electron- 
hole and gas plasma, semiconductor and superconductor 
structures, systems with uniformly generated combus- 
tion material, autocatalytic chemical reactions, models 
of population dynamics P-pJ|. Recently, an intriguing 
phenomenon was observed in the chemical experiments 
with Turing patterns |l2[ , combustion system with cel- 
lular flames and the catalytic reaction on the sur- 
face |Q. In these experiments the patterns that were 
observed consisted of many disk-shaped domains whose 
radii oscillated in time. The oscillations of different do- 
mains either synchronized or exhibited complex spatio- 
temporal behavior. 

Many nonequilibrium systems in which patterns can 
form, including the systems mentioned above, are de- 
scribed by the reaction-diffusion systems of the activator- 
inhibitor type, the simplest of which is a pair of reaction- 
diffusion equations [p|-p"l| 



re 



do 



I 2 AO - q(9, ry, A), 



(1) 



(2) 



where is the activator, r\ is the inhibitor, I and L are the 
characteristic length scales and Tg and are the char- 
acteristic time scales of the activator and the inhibitor, 
respectively; q and Q are certain non-linear functions, 
and A is the system's excitation level. For example, in 



the system with the uniformly generated combustion ma- 
terial the activator is the temperature of the gas mixture, 
the inhibitor is the density of the fuel, A is proportional 
to the rate of the fuel supply, and the non-linear func- 
tions q and Qcontain the dissipation, supply, and reac- 
tion terms |||9) . The length scales / and L are related to 
the diffusion coefficients of the activator and the inhibitor 
Dg and D v , respectively: I — \JUgTg and L = ^JD^t^. 

Kerner and Osipov showed that the properties of the 
patterns forming in the systems described by Eqs. (Q) 
and (|^) are determined mainly by the parameters e = l/L 
and a — rg/r^ and the form of the nullcline of Eq. ( |l|) . 
For many systems this nullcline is N-shaped (Fig. Tlh 
^[-|l0|. In such N systems static domain patterns form 
when e< 1 and a 3> e (KN systems), traveling waves 
(autowaves) at a <C 1 and a < e 2 (SIN systems), and 
all sorts of dynamic patterns at e <C I and e 2 < a < e 



(KSIN systems) ||-|§(]JJ]J. As a result of the instabil- 
ity of the homogeneous state periodic and more complex 
patterns form in these system, whereas when the homoge- 
neous state of the system is stable one can excite solitary 
patterns — autosolitons (AS) — by a sufficiently strong 
localized external stimulus. In KN systems these pat- 
terns are collections of static domains with sharp walls 
(interfaces) whose width is of order I J|,|].[l5]jl(]] . These 
domain patterns may undergo different kinds of instabil- 
ities leading to the formation of complex dynamic pat- 
terns when a becomes sufficiently small. The simplest 
example of such destabilization is the transformation of 
a static AS to a pulsating (breathing) AS. This effect 
was discovered by Koga and Kuramoto in an axiomatic 
reaction-diffusion model JItJ and subsequently studied 
by many authors, both for one-dimensional and higher- 
dimensional radially-symmetric AS [^]-[l0 15 T^Jf^]. Pul- 
sating AS were directly observed in semiconductors p0[ , 
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composite superconductors pH , autocatalytic reaction 
p2| , and combustion experiments (2^] . 

The situation becomes much more complicated when 
instead of a single AS the pattern consists of many inter- 
acting domains. Several attempts to approach this prob- 
lem were made for one-dimensional N systems. Kerner 
and Osipov showed that the stationary periodic domain 
patterns (strata) may undergo instability and transform 
into a breathing pattern 10|. Ohta et al. were able 
to obtain the linearized equation of motion for the pe- 
riodic patterns in the piecewise-linear reaction-diffusion 
model and showed that the introduction of a simple non- 
linearity to these equations leads to the synchronization 
of the domain oscillations ]24| , |2"5| ] . Still, the question of 
the effect of the interaction of domain patterns under- 
going breathing motion, especially in higher-dimensional 
systems, remains largely unresolved. 

It is easy to show that in the case t « 1 the wave 
length of the fluctuation with respect to which the Turing 
instability of the homogeneous state of the system (|l|) 
and (|) is realized is A - (LL) 1 / 2 |@. For this reason 
the periodic Turing structures that form in the system 
have the period C p ~ e x / 2 L <C L. Moreover, in the 
higher-dimensional systems with t«l any static pattern 
will consist of the domains whose characteristic size is 
of order e^ 3 L < L Jl|,|l|,||. This means that in a 
pattern consisting of many domains one domain interacts 
with a large number of other domains at the same time. 
This fact should significantly reduce the complexity of 
the interactions between different domains. 

In this paper we will study the periodic one- 
dimensional and two-dimensional hexagonal patterns of 
small period (C p -C L) undergoing the oscillatory in- 
stability in an arbitrary Kf2N system. Using the inter- 
facial dynamics approach, we will derive the universal 
nonlinear equations describing the pulsations of the peri- 
odic multidomain patterns in arbitrary dimensions. We 
will analyze the conditions for the synchronization and 
breakdown of the pulsations, the effects of the disorder, 
and study possible complex spatio-temporal behaviors. 
Finally, we will compare the effects studied by us with 
relevant experiments. 

Our paper is organized as follows. In Sec. II we reduce 
Eqs. ([!]) and (^) to the problem of interfacial dynamics 
in the limit e — > and show the way of treating this 
problem in the case C p <C 1. In Sec. Ill we apply this 
method to the one-dimensional periodic strata. In Sec. 
IV we consider two-dimensional hexagonal patterns and 
also briefly discuss the case of the three-dimensional pat- 
terns. In Sec. V we analyze the domain oscillations in a 
disordered pattern with a random distribution of domain 
sizes. In Sec. VI we present the results of the numeri- 
cal simulations for a concrete model, and, finally, in Sec. 
VII we discuss the relevancy of the obtained results to 
the experiments and draw conclusions. 



II. INTERFACIAL DYNAMICS PROBLEM 

Let us measure the length and time in the units of L 
and t v , respectively. Then Eqs. (]l|) and (||) become 



a—=e 2 A6-q(9,r,,A), 



^ = Ar,-Q(9, V ,A). 



(3) 



(4) 



The boundary conditions for these equations may be neu- 
tral or periodic. From the mathematical point of view the 
fact that 9 is the activator and r\ is the inhibitor means 
that for some values of 9 and r\ we have q' e < 0, and that 
for all 9 and r) 



Q; > 0, q' v Q' e < 0, 



(5) 



and the derivatives Q' vl Q' ei and q' v do not change sign. 

The dynamics of the domain patterns forming in KN 
and Kf2N systems can be reduced to the interfacial dy- 
namics problem in which the dynamics of the pattern 
interfaces is coupled to the bulk field Jl8|,^6| . Far enough 
from the domain interfaces (at distances much greater 
than e) the bulk field must satisfy the equation of smooth 
distributions (outer solution) 



^=Ar l -Q(9(r ] ),r,,A), 



(6) 



in which 9 and r\ are related by the equation of local 
coupling 



q(9, V ,A) = 0, 



(7) 



that is, far from the domain interfaces 9 and r\ lie on the 
nullcline of Eq. (||). This relation is multivalued (see 
Fig. [I]), so one has to choose the branch with 9 < 9 in 
the "cold" regions (the domains of low values of 9) and 
9 > 9' in the "hot" regions (the domains of high values 
of 9). The dynamics of the interface is governed by the 
following equation 



Of 
St 



2 a 1 K(r) + v(r)(f)), 



(8) 



where f is a point on the interface, n is the normal vector 
to the interface pointing into the cold region, K(r) is the 
curvature at the point r of the interface, and u (17(f)) is 
the velocity of the interface in the absence of the cur- 
vature as a function of the value of the bulk field 77(f) 
at the interface. The function v{rf) is a solution to the 
nonlinear eigenvalue problem and is in general some com- 
plicated nonlinear function of 77 |T^ , ^| . 

The interfacial dynamics problem represented by Eqs. 
(||) and (|J) is a highly nonlinear problem. However, the 
situation becomes simpler if the pattern consists of the 
alternating hot and cold regions whose characteristic size 
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is much smaller than 1. Kerner and Osipov showed that 
in this case the value of 77 is close to a constant in the 
entire pattern The reason is that because of its 

high diffusivity the inhibitor cannot react well on these 
variations of the activator whose characteristic length is 
much smaller than the characteristic length scale of the 
inhibitor. This allows us to linearize Eq. (^) around 
rj = r) s , where r/ s is the constant value of the inhibitor. 
Introducing 



(9) 



we get 
di) 

~di = 
where 



A77 - cifj - (c 3 - ci)fjl(x) ~ Q(6si,r)s) - al(x), 

(10) 



Cl,3 = Q v (0 s l,3,Vs) '- m r , (11) 



a = Q(0 S 3,Vs) - Q{0 s i,t] s ), 
the values of # s i,3 satisfy 

q(0 s i,3,Vs) = 0, 



(12) 



(13) 



8 s i and 6 S 3 being the minimal and the maximal roots 
of Eq. (U|), respectively; I(x) is the indicator function 
which is equal to 1 if x is in the hot region and zero in 
the cold region. Notice that, according to Eqs. (g), the 
constants c\ and C3 are positive. 

Obviously rj s should be equal to the value of 77 at which 
v(r}) = in order for the pattern to be stationary. The 
value of 77 s must therefore satisfy [p|,|i"o|,p6[ 



q(9,r) s )d0 = 0. 



(14) 



For small fj the function 11(77) may be linearized around 
T7 S , so we obtain Q 



where 



and 



eBfj 
otZa ' 



B = -a q' v (e,vs)de, 



(15) 



(16) 



V-2U e cL6, U e = - q(0, Vs )d8. (17) 



Z = 



The constants B and Z are of order 1. Notice that, ac- 
cording to Eq. (||), the value of B is positive. 



Although Eqs. (|8j) and ( |10D with (|15|) are simpler than 
the original interfacial dynamics equations, they are still 
difficult to deal with. It appears, however, that these 
equations can be further simplified for treating multido- 
main patterns by introducing some averaged variables. 
We will outline this procedure in this section and demon- 
strate its application to the periodic patterns in the sub- 
sequent sections. 

Let us introduce the number of domains per unit vol- 
ume n and the radius of the domain p (in one dimension 
p is the half- width of a domain) . If we now average Eq. 
( |Io| ) over the volume of size s such that p < s <C 1, we 
will obtain the "coarse-grained" equation for the average 
value of the inhibitor 77 = (fj) 



^ = A?7 - ci7j(l + c 2 np d n d ) 



,1,77s) - anp d VLdi 
(18) 



where fid is the volume of the d-dimensional unit sphere 
(ill = 2), and 



(-2 



= (C3 - Cl)/Cl 



(19) 



measures the asymmetry between the hot and the cold 
domains. According to Eq. (|5|), the value of c 2 > — 1. 
Besides 77, there is a local contribution to 77 due to the 
variation of 77 on the length scales of p 1. Let us 
introduce 



77 = 77 - 77 



(20) 



Since 77 varies on the short length scales, the terms pro- 
portional to 77 in the right-hand side of Eq. (|Io| ) are small 
compared to the Laplacian as well. Also, when the char- 
acteristic time scale of variation of 77 is much greater than 
p 2 , as is the case in all interesting situations (see the fol- 
lowing sections), the time derivative in Eq. (|10|) is small 
compared to the Laplacian. Subtracting Eq. (|18| ) from 
Eq. (|l(]) and neglecting all these terms, we obtain 



A77 = a{I(x) - (I(x))}, (fj) = 0. 



(21) 



From this equation one can obtain the value of 77 = 7% in 
the wall of an individual domain, so from Eqs. (^) and 
( |l5| ) follows that the equation for the radius is 



dp 

di 



e 2 (d-l) eB{fj + fj s ) 



OLp 



otZa 



(22) 



The variables 77, 7%, and p may now be considered as 
space- and time-dependent on the length scales much 
greater than s. Thus, we have a closed set of equations 
for these coarse-grained variables, so the number of the 
relevant dynamical variables in the problem is consider- 
ably reduced. 
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III. ONE-DIMENSIONAL PERIODIC STRATA OF 
SMALL PERIOD 

Let us consider the one-dimensional periodic strata of 
the period C p <C 1 (Fig. ||). In this case n = C~ l , so Eq. 
(fL8h becomes 



§ = - C!fj(l - 2c 2 £;V) - ^L-\ P - p ), (23) 
where we introduced 

p o = -C p Q(0 a i,ri s )/2a. (24) 



Equation (|21j) for a one-dimensional pattern with the pe- 

C( 

d 2 fj 



riod C p becomes 



i ,-al(x)-2a£- 1 p, (fj) = 



(25) 



with neutral boundary conditions at x = ±£ p /2 (because 
of the translational invariance, we may choose the cen- 
ter of the domain to be at i = 0). A straightforward 
calculation gives 



Vs = a 



6 ' 3£„ 



Rescaling the variables r\ and p 

fj' = f)/aC 2 p , p = p/C p , 
and introducing the quantities 



(26) 



(27) 



n = aZ/eBC p , uj = ^2eB/aZ£ p . (28) 

we write Eqs. ( |22] ) and (^3|) (dropping the primes) as 
dfj d 2 fj 



dt dx 2 



Cl fj(l + 2c2p)-^Ti(p-po), (29) 



dp _ p 2 4p 3 



(30) 



The parameters of the original reaction-diffusion system 
[Eqs. (||) and (0)] enter these equations only through 
constants Tb and c 2 (the value of ci can be absorbed 
in the definitions of ri and ojq), and the excitation level 
enters mainly through p~Q. Note that since, according to 
its definition, the value of 2p cannot exceed £„, we have 
< p < \. Also note that Eqs. (§|) and @ with the 
spatial derivative equal to zero describe the oscillations 
of a single domain in the system of size £ p < 1. 

Equations of the type of Eqs. ( |29| ) and ( J30| ) have been 
extensively studied in the context of oscillatory chemical 
reactions P-p|,p7| . It was shown that these systems host a 
great richness of dynamic patterns, such as self-sustained 
uniform oscillations, traveling waves, target patterns and 
spirals (in higher dimensions), and spatio-temporal chaos 



(turbulence). This is also true for Eqs. (|28|) and (|3£ 
However, one should keep in mind that these equations 
describe the dynamics of the already existing multido- 
main pattern, in other words, the dynamical patterns 
described by Eqs. (^9|) and ( |30| ) will be seen "on top" of 
the stationary multidomain (Turing) patterns. This is an 
important distinction from the oscillatory systems, such 
as oscillatory chemical reactions, in which the dynamic 
patterns appear on top of the stationary homogeneous 
state. 

We are interested in the collective oscillations of the 
domains in the pattern. According to Eqs. ( p9| ) and 
(p0|), the characteristic frequency of these oscillations is 
equal to too. In view of Eq. ([28]), the period of oscillations 
is the shortest time scale in the problem if 



(31) 



Recall that in deriving Eq. (pil| ) we neglected the time 
derivative of fj. This is justified if uj a <C £p 2 , what is 
equivalent toa> eCp~, so the first condition in E q. (|3l| ) 
is in fact a necessary condition for Eqs. ( p9| ) and (3^) to 
be valid. 

Equations (|29| ) and ( |30| ) have a trivial solution p = 
const, fj = const which corresponds to the stationary 
multidomain pattern (strata). For a satisfying the con- 
dition in Eq. @ we have luq 1, so the value of p 
for the stationary pattern is close to po. Linear stability 
analysis of Eqs. (|29| ) and (|30| ) shows that at t\ > ri c , 
where 



Tu 



12p - 2Ap 2 - 1 
6c 



and 



c = ci(l + 2c 2 po) > 0, 



(32) 



(33) 



this solution is stable. At t\ — t\ c it looses stability with 
respect to the uniform oscillations with the frequency 
u = u>q and the wave vector k = 0. The uniformly os- 
cillating solution corresponds to the pattern in which all 
domains are oscillating in phase with the frequency luq. 
This is a supercritical Hopf bifurcation (see below). In 
view of Eq. fl28y and d32h, the value of n c ~ 1, so in 



terms of the original variables we have a c ~ eC p and 
Also, according to Eq. (|32|), the value of 



i~,p . 



Tic > only if 

1 

4 



1 

7! 



< Po < 



1 

71 



(34) 



i.e., 0.11 < po < 0.39. This is a completely universal 
result for all one-dimensional strata of small period in 
Kf2N systems. It is easy to see from Eq. ( p2[ ) that if the 
condition of Eq. (|34|) is satisfied, for any acceptable value 
of c 2 there exists T\ c = T^ ax , such that the stationary 
pattern is stable for n > r™ ax for any value of pq. If 
the condition in Eq. (|34|) is not satisfied, the pattern will 
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still become unstable when a ~ e 2 1| [n],[l^]. This will 
happen at the values of a which are much smaller than 
in the former case. Indeed, according to Eq. (|3^), the 
oscillatory instability occurs at a > e 2 since C p is always 



much greater than e [|B|,[10|. 

To study the destabilization in more detail we will use 
the fact that the period of the oscillations Uq 1 is much 
smaller than the characteristic relaxation time n , so one 
can a pply the method of Bogoliubov and Mitropolsky to 
Eqs. © and @ §|. To do this, we will use Eq. @ 
to eliminate fj from Eq. (E97) and write 



po + Rcos(Lu t + 6), 



(35) 



where R and O are slowly varying functions of t and x. 
Substituting Eq. (|35|) into Eq. (p9|), multiplying it by 
sin(wo^ + 0) and cos(woi + 0) and integrating over the 
period, we will get, respectively, 

dR cR, . R 3 ld 2 R RfdQ\ 2 ,„„, 

-m=^ inc - Tl) -2^ + 2d^- 2 Ite) ' (36) 



and 



i? 9x 9a; 



(37) 



In deriving Eqs 
terms. 

Equations 



and 



we kept only the leading 



(36) and 



ft are equivalent to the equa- 
tions for the amplitude and phase for the well-known 
Newell- Whitehead equation |29|| . An important class of 
solutions of this equation is R = Rq and O = kx, where 



Ro = y c(t\ c - n) - rife 2 



(38) 



and k < \/c(j\ c — t\)Jt\. It represent the nonlinear dis- 
persionless traveling waves of the phase of the domain 
oscillations with the wave vector k which are stable if 
k < \J c(t\ c — t\)/Zt\. A particular solution with k = 
corresponds to the synchronous oscillations of the do- 
mains. In fact, if the boundary conditions for the origi- 
nal reaction-diffusion system (^) and (0) are neutral, as 
in most of the real situations, this is the only admissi- 
ble plane wave solution. Since it is stable, a variety of 
initial conditions will evolve into it, causing the oscilla- 
tions of different domains to synchronize. This means 
that the synchronization of oscillations of different do- 
mains is the major scenario of the development of the 
small-period patterns in one dimension. Equations j3t^ ) 
and ( |37[) also admit a solution in the form of a domain 
boundary (kink) , upon going through which the phase O 
changes by it. These are the only stable stationary solu- 
tions of Eqs. ( |36| ) and (|37|), so any initial condition will 
in general produce a collection of regions in which the 
domains oscillate in phase. In the regions next to each 
other the domains will oscillate in antiphase. At t\ < t\ c 
Eqs. ( p6| ) and ( |37| ) also have the solutions in the form of 
the wave in which the stable state R = Rq with k = 



and = const invades the unstable state R = 0. This 
wave corresponds to the onset of the synchronous oscil- 
lations in the unstable stationary strata from a localized 
initial condition. The propagation velocity of this wave 

is v* = y/c(r lc - n)/ n @. 

It is clear that small local inhomogeneities, both in- 
trinsic and extrinsic, may play the role of the organiz- 
ing centers for the one-dimensional analog of the target 
patterns: waves traveling away from the inhomogeneity 
in both directions. An example of an intrinsic inhomo- 
geneity here is a non-uniform distribution of the domains 
along the x-axis. In this case instead of the constant den- 
sity of domains we have a smooth function n(x), which 
can also vary in time and can be considered as an ex- 
tra degree of freedom in the problem. The density must 
satisfy the continuity equations 



where 



m + dx {nY) = °> 



p drj 
t 2 dx 



(39) 



(40) 



is the velocity of the domain as a whole due to the gra- 
dient of fj [see Eq. (|l5|)1, and 



T 2 = aZ/eBCl 



(41) 



is the characteristic time scale of the density variation. 
For Ti ~ 1 we have r 2 ~ £~ 2 , so the density relaxes on 
much longer time scale and any non-uniform distribution 
of n is frozen on the time scale of the relaxation of the 
amplitude and phase. This means that by setting a non- 
uniform density of domains as the initial condition one 
can produce complex spatio-temporal patterns of the do- 
main oscillations which will not synchronize for a very 
long time. 

From the definition of R follows that it lies in the in- 
terval from po to h — po- If R exceeds one of these values, 
the domains will either merge or collapse in the course 
of the oscillations. On the other hand, according to Eq. 
(p8|), the value of Rq may exceed either po or | — po 
for some values of t%. This can be seen from Fig. |^, in 
which the bottom solid line corresponds to the values of 
t\ at which this happen. Figure || thus shows that the 
synchronization of the domain oscillations occurs only in 
the relatively narrow region of the system's parameters. 
The merging or collapse of the domains will result in the 
breakdown of the collective domain oscillations and col- 
lapse of the strata. Indeed, suppose that at some moment 
every two domains merged into one. This may be viewed 
as doubling of the pattern's period L p . According to Eq. 
( p8| ) , this will result in the decrease of T\ and further in- 
crease of Rq. Therefore, the process of domain merging 
will have an avalanche character and lead to a complete 
destruction of the strata. The way the breakdown of the 
strata will follow will depend on whether or not the ho- 
mogeneous state of the systems is stable with respect to 
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the uniform self-oscillations of the homogeneous state at 
this value of A. If the homogeneous state of the sys- 
tem is unstable, the breakdown will result in the onset of 
the uniform relaxation self-oscillations (the latter follows 
from the fact that since t\ c < 1 and e « 1 we have a«l 
and, therefore, relaxation oscillations §-0). If the ho- 
mogeneous state of the system is stable, the breakdown 
will result in the formation of traveling domain patterns 
(AS) or in the collapse of the strata into the homogeneous 
state. 

The scenarios of the evolution of one-dimensional 
strata of small period discussed in the previous paragraph 
are precisely what we see in the numerical simulations of 
a concrete model. We also find a detailed quantitative 
agreement between the limit cycle oscillations of a single 
domain determined by Eqs. (|2^) and ( p0| ) and the re- 
sults of the numerical simulations of a reaction-diffusion 
model for small enough e and C p . 

Before concluding this section, let us discuss another 
dynamic behavior of strata in one dimension. As was 
recently shown by Osipov, for a < e periodic patterns 
may undergo the instability leadingto the formation of 
the traveling strata (a wave train) pl[ | . Osipov obtained 
the general criterion for this instability. In the case of 
one-dimensional strata of small period it is possible to 
calculate the threshold value of a, since the spatial vari- 
ation of 77 is given by i)(x). Substituting it into the cri- 
terion (Eq. (4.5) of Ref. [Q), we obtain that the trans- 
formation of the stationary strata to traveling occurs at 
ctT ~ eCp One can see that this instability occurs for 
much smaller values of a than the oscillatory instability 
for which a c ~ eC p . Therefore, the instability leading 
to the formation of the traveling pattern can never be 
reached for the one-dimensional strata of small period. 

Throughout the analysis presented above we used the 
fact that the period of the pattern is small compared to 
the length of the variation of the inhibitor and used C p 
as the small parameter in the expansions. The smallness 
of the parameter £ p greatly simplified the treatment of 
the domain interactions and lead to the universal results 
which are practically independent of the concrete non- 
linearities of the system. This is a general property of 
the KN and KON systems [|5[^6|j2|] . Unfortunately, no 
such treatment is possible for the one-dimensional peri- 
odic patterns whose period is of order 1. In the latter 
case one is faced with the problem of treating the dy- 
namics of each individual domain separately and taking 
into account the complex memory effects associated with 
the domain interactions, which cannot be solved in gen- 
eral. So, the problem of the interaction and dynamics of 
the one-dimensional patterns consisting of the domains 
whose size if of order 1 remains open. This problem, 
however, does not exist in higher dimensions, if e is small 
enough. In this case the domain sizes are necessarily 
small since the domains whose size is greater than e 1 / 3 
are always unstable with respect to the transverse insta- 
bility of their walls and split into smaller domains fl5|.p6| . 



IV. HIGHER-DIMENSIONAL MULTIDOMAIN 
PATTERNS 

Let us now study periodic multidomain patterns in 
higher dimensions. For definiteness we will consider a 
two-dimensional hexagonal pattern of disk-shaped do- 
mains, the most frequent Turing pattern in chemical ex- 
periments [^|. As was shown by Muratov and Osipov, in 
two dimensions a stable stationary multidomain pattern 
must consist of the domains whose size is of order e 1 / 3 
Jl5| , [l6| . Since the period of the pattern is of the same or- 
der as the domain size, it can be used as a natural small 
parameter, so all the ideas of the analysis of the previous 
section should apply to the two-dimensional multidomain 
patterns. 

Let us proceed with the derivation of the equations 
describing the collective oscillations of the domains in 
the pattern. For a hexagonal pattern with the period C p 
the density of the domains is n — 2/^C p , so Eq. JlS| ) 
becomes 

d 4 = Afj - cfj ( 1 + - ~ M , (42) 



where we introduced 



Pi 
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(43) 



Because of the anisotropy of the hexagonal lattice the 
shape of the domains will in general slightly deviate from 
the disk shape. It is clear that this deviation is small 
and is smaller when the radius of the domains is smaller. 
So, for the sake of simplicity and without significantly 
affecting the results of the analysis, we will ignore the 
effects of the anisotropy and consider the domains to be 
ideally circular. Also, instead of considering the problem 
for 77 in the hexagonal cell, we will solve it for a circular 
domain of radius C p /2. Thus, we have 



d 77 
dr 2 



1 dfj 
r dr 



= al(r) 



Aap 2 
p 



(v) = 0, 



(44) 



where r is the radial coordinate and the boundary con- 
ditions are neutral at r = and r = C p /2. Solving this 
equation for 77 and calculating f/ s = t)(p), we get 



rj s = a 



3p 2 3p 4 



so, Eq. rt2 



dp 
dt 
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becomes 
eB 



aZa 



2C 2 



8 



P! ln 2p 
2 £ T) 



2CI 



2 L,nQ 



(45) 



(46) 



Rescaling the variables according to Eq. ( p7| ) and drop- 
ping the primes, from Eqs. (E2n and (|46]) we get 
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(47) 



(48) 



to the fluctuation with u> = loq and k = when ri < t± c , 
where 



4c + 5pg-24pg + 4pgln2po 



(52) 



Going back to the original variables, we see that in the 
case of the two-dimensional multidomain pattern the 
destabilization occurs at 



e = eZ/BCl- 



(49) 



,4/3 



(53) 



As before, < po < i. The parameter e does not ap- 
pear in the one-dimensional case. It characterizes the 
stability of the domains with respect to the transverse 
perturbations. Since the curvature radius of the stable 
domains must be of order e 1 / 3 [fi5| , [l6| , the value of e is 
of order 1 for the domains whose size is comparable with 



r 



A/3 



Notice that Eqs. (|7| and (fcfwithout the 



space dependence and with the coefficient ir/ y3 replaced 
by 2 in Eq. (47) will describe the oscillations of a single 
domain in a circular region of radius C p <C 1 . This situ- 
ation was realized in the experiment by Haim et al. on 
the ferrocyanide-iodate-sulfite (FIS) reaction p2fl . 

When the condition of Eq. (^l|) is satisfied, Eqs. ( p7| ) 
and (Q) have two time scales, as in the one-dimensional 
case: the short time scale uJq 1 for the oscillations and the 
long time scale n for the relaxation of their amplitude 
and phase. However, in contrast to the one-dimensional 
case, the analysis is complicated by the fact that now the 
oscillations themselves are nonlinear. To proceed further, 
to express fj in terms of p and sub- 
47|), taking into account that uj ^ 1- 



we will use Eq. (|4£ 
stitute it into Eq 
We obtain 

d 2 p 7T 2 2 



Po) 



Tl 



A ) — , 

at' 



(50) 



where 



/(p) = c i T i 



4^-^+6p 3 -p\n2p. 



(51) 



The left-hand side of Eq. ( p0| ) is an equation of mo- 
tion for a particle of unit mass in the potential V = 
^8- ^ — pQp^J . For < p < \ it describes finite motion 

between p = p m i n and p — p max with the characteristic 
frequency uq. The right-hand side of Eq. ( ^o|) is a weak 
nonlinear friction force which also contains the diffusion 
term that couples the oscillators in a spatially distributed 
system. 

As in the case of the one-dimensional patterns, the 
equilibrium solution p — p~o of Eq. ( pp[ ) corresponds 
to the stationary multidomain pattern. Equation j50| ) 
shows that the stationary multidomain pattern becomes 
unstable (the friction / becomes negative) with respect 



Note that the value of a c is smaller than a u ~ e at which 
a localized domain (AS) is destabilized with respect to 
pulsations ||j||l5|. 

Figure [| shows the possible forms of the dependences 
T~ic{po) obtained from Eq. (|5^ ) for different values of e. 
One can see from Eq. (52) that for e > 0.031 for any value 
of po there exists a value of t± at which the instability 
occurs. When 0.00045 < e < 0.031 the instability is 
realized only for p < p a3 with p a3 < |, whereas for 
e < 0.00045 the instability is realized for p < p al and 
Pai < Po < Pa3 (see Fig. |). When e > e* (for example, 
e* = 0.0034 for C2 = 0) for a given n the instability is 
realized for a single value of po, whereas for e < e* there 
are values of t\ for which the instability is realized for 
three values of po- 

Let us take a closer look at Eq. (|50|). The behavior 
of the oscillations is determined by the sign of the fric- 
tion term / in Eq. ([sb]). Three possible forms of this 
term as a function of p are shown in Fig. [|. The os- 
cillation amplitude decreases if p remains in the domain 
where / > 0, or increases if p is in the domain where 
/ < (indicated by arrows in Fig. |^) . In the situation of 
Fig. ||(c) all values of p correspond to the amplification 
region, so the amplitude of any oscillations will increase 
until p m ax becomes equal to h or until p m in = 0. In the 
first case the neighboring domains will merge while in 
the second the domain will collapse. Both these effects 
will cause the restructuring of the pattern and significant 
changes in the collective domain oscillations. In the sit- 
uation of Fig. ||(b) the oscillations arc amplified when 
p < p C 3 and attenuated when p > p c ^. The stationary 
multidomain pattern will be stable when po > p c % or un- 
stable when po < p C 3- Depending on the signs of the 
first and the second derivatives of f(p) the bifurcation at 
Po = Pc3 may be both supercritical and subcritical (this 
is a Hopf bifurcation). It is clear, however, that deeper 
into the unstable regime the oscillations will be amplified 
until they collapse when p m i n — even in the case of the 
supercritical bifurcation. 

The most interesting situation is shown in Fig. ^|(a). 
There the oscillations are amplified when p < p c \ or 
Pc2 < P < Pc3 and attenuated for the rest of the val- 
ues of p, so the pattern is stable only for p c \ < po < p C 2 
or po > Pci- The bifurcation at p c ^ is supercritical, while 
the bifurcation at po = p c i is subcritical. This can be 
easily seen from the following argument. By varying the 
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value of T\ one can make the minimum of / between p C 2 
and pc$ arbitrarily shallow, thus controlling the ampli- 
fication of the oscillations in the region p C 2 < p < p c z- 
Since in this case the small region of amplification is sur- 
rounded by a finite region of dissipation, the bifurcation 
at p C 2 and p C 3 will be supercritical. Since the first and 
the second derivatives of / at p c i are opposite to the one 
at p C 2, the bifurcation at p c \ is subcritical. If the value 
of Ti is decreased, the minimum of / between p C 2 and 
Pci gets deeper, and at some value of t\ the bifurcation 
at p C 2 becomes subcritical. Upon further decreasing the 
value of T\ the limit cycle associated with the subcritical 
bifurcation at p C 2 will become unstable and lead to the 
collapse of the oscillations. Thus, the most stable limit 
cycle oscillations should be expected when p C 2 < po < p c z 
and the minimum of / between p C 2 and p C 3 is sufficiently 
shallow. The analysis of the bifurcation types for dif- 
ferent values of e and po is summarized in Fig. ^ (for 
simplicity we used C2 = 0). 

Up to now we discussed the oscillations of a single do- 
main in the multidomain pattern. What we found is 
that in certain situations the synchronous oscillations of 
all domains in the pattern can be stable. This natu- 
rally implies that the synchronously pulsating multido- 
main pattern is an attractor and the synchronization of 
the domain oscillations is therefore one of the major sce- 
narios of the development of the stationary multidomain 
patterns. However, the situation is reacher in two dimen- 
sions because of the possibility of other classes of stable 
solutions, such as target patterns and spiral waves. In 
principle, it is possible to use the fact that u^ 1 <C T\ 
and obtain the equations for the amplitude and phase by 
averaging over the trajectories of the frictionless oscilla- 
tors in this case as well. This calculation however, can- 
not be carried out in the analytical form. On the other 
hand, in a wide region of the parameters the bifurcation 
of the stationary pattern is supercritical, so in its vicinity 
the dynamics of the pattern is described by the complex 
Landau-Ginzburg equation. In the appropriately scaled 
time and space variables this equation in the considered 
case has the form 

^ = AW + W - (1 + ib)\W\ 2 W, (54) 
ot 

where W is appropriately normalized complex amplitude 
and 



b = 



4 • 2V2 . 33/4^1/2 ,5/2 



Pq UpTic 



^c lC2 T lc p 4 - 56e + 216p§ + p* + 4p£ In 2p Q 



(55) 



The complex Landau-Ginzburg equation has been the 
subject of intensive studies for the last two decades (see, 
for example, Refs. P, |27| , p2| , ^3[ ) . Specifically, for the case 
of Eq. (Q) Hagan obtained the solution in the form 
of a steadily rotating spiral wave j34j ]. In our case this 
spiral wave will be seen on top of the stationary multido- 
main pattern. Precisely this phenomenon was observed 
by Boissande, Dulos, and De Kepper in the experiments 



with the Turing patterns in the chlorite-iodide-malonic 
acid (CIMA) reaction |l2"||. 

Kuramoto and Koga showed that the spiral waves in 
Eq. (|54| ) become unstable at sufficiently large b 
They showed that the spiral breakup leads to the for- 
mation of a chaotic spatio-temporal pattern — spiral 
turbulence. This will be the case in the situation we 
study. Indeed, throughout the analysis we assumed that 
uqTi 1, what means that b ^> 1, so in general the 
spiral wave solution is unstable. Note, however, that the 
value of b decreases as C p increases, so when C p (and e) is 
not very small, the value of b may become small enough, 
so that the spiral wave solution is stable. This can also 
be achieved by making the coupling constant B smaller. 

Since generally b 1, plane wave solutions will be 
stable only in a narrow range of the wave vectors around 
k = |27]] . This is because, in contrast to the one- 
dimensional case, the dispersion of the plane waves is 
high. However, for Eq. ([34]) the Benjamin-Feir insta- 
bility is not realized, so the synchronous oscillations of 
the domains are always stable in the spatially distributed 
system close to the bifurcation point. 

The stationary multidomain patterns are never per- 
fect. The defects of the pattern will work as initiators 
for the spiral waves, which will in turn break up and pro- 
duce more spiral waves invading the entire pattern, so 
the formation of the chaotically oscillating multidomain 
pattern is the more likely scenario for the development 
of the unstable stationary multidomain pattern in two 
dimensions. 

The imperfections of the multidomain pattern may also 
work as the guiding centers generating target waves of 
the oscillation phase j27|. An argument similar to the 
one applied to the one-dimensional strata may be ap- 
plied here. The target waves may be initiated by the 
smooth spatial variations of the density of domains. The 
equation for the density will be given by Eqs. (39) and 
(pio|), if one replaces d/dx by V, so the time scale of the 
density variation will once again be much longer than n , 
and, therefore, the density can be considered frozen on 
the latter time scale. However, since the dispersion of 
the waves in the considered situation is high, the target 
patterns are likely to be unstable, what should also lead 
to the stochastization of the domain oscillations. 

It is clear that qualitatively the scenarios discussed 
above are also realized when t% is not very close to t\ c . 
Of course, for arbitrary t\ one should also take into ac- 
count the possibility of the collapse or merging of the 
domains in the course of the oscillations. To determine 
the values of t± at which the collapse or merging of the 
domains occur for a given value of po one can consider the 
oscillator in Eq. ( |50|) without the space dependence. The 
collapse or merging of the domains occurs when p m i n = 
or p max = |) respectively. The values of p mm and p max 
can be obtained from the condition of energy balance 
for the steady oscillations. Recalling that the friction is 
weak, from the elementary mechanics we get 



8 



where is the energy of the oscillator corresponding to 
the values of p m in and p ma x- We solved this equation 
numerically in the case C2 = 0. Figure |?] shows a typi- 
cal dependence of the amplitude of the oscillations on T\ 
for a given value of e in the case of the supercritical bi- 
furcation (pmin and Pmax correspond to the left and the 
right portions of the solid curve, respectively). From this 
figure it can be seen that the limit cycle oscillations are 
stable only if T\d < t± < t\ c . If t\ < T\d, the amplitude 
of the oscillations will increase until the domain collapses 
or merges with the neighbors. Several possibilities exist 
for the values of p min and p max - The value of T\d may 
happen to be negative, in this case the stability of the 
limit cycle oscillations will depend on whether the value 
of Pmax becomes greater than ^ at some value of T\ . The 
oscillations will be stable for t\ greater than this value if 
this is the case, or stable for all values of t\ in the op- 
posite case. The value of p max may also happen to be 
greater than | . In this case the region of the stability of 
the domain oscillations will also be narrower. Notice that 
Pmin i s always greater than zero because of the singular 
character of f(p) at p = 0. 

In the case of subcritical bifurcation the limit cycle os- 
cillations are unstable for any values of ri, so collective 
oscillations of the domains in the multidomain pattern 
always break down in this case. The most likely sce- 
nario here is that some of the domains will shrink and 
collapse whereas some will grow, what can be effectively 
thought of as an increase of the pattern's period. As a 
result, the value of e, which mainly determines the type 
of the bifurcation, will get smaller, so in the new pat- 
tern of greater period the bifurcation may happen to be 
supercritical. This means that the pattern may natu- 
rally evolve into a state in which the collective domain 
oscillations are stable from more or less arbitrary initial 
conditions, provided that the initial condition consists of 
the alternating hot and cold domains of size less or of or- 
der e 1 / 3 . It is also possible that the pattern that forms in 
this process is stable with respect to the domain oscilla- 
tions. Yet another possibility is that in this process some 
of the domains will merge, what will result in the for- 
mations of an irregular pattern, which also may exhibit 
chaotic dynamics. The merging process may also lead to 
the destruction of the multidomain pattern. We will get 
back to these points when we will discuss the results of 
the numerical simulations of a concrete model. 

In the analysis presented in this section we did not 
specify precisely the range of the values of e. This range 
is determined by the stability of the stationary multido- 
main pattern with respect to the activator repumping 
and the transverse distortions of the walls of the do- 
mains. When e is greater than some value for a given 
Po, the pattern becomes unstable with respect to the ac- 
tivator repumping effect leading to the doubling of the 
pattern's period, whereas when e is smaller than some 



critical value at a given po, the domains destabilize with 
respect to the transverse instability of their walls leading 
to the radially-nonsymmetric distortions of the domains 
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Both these instabilities occur when e ~ 1 in 
the limit e — > regardless of the value of a. It is possible 
to show [^6| that for the stable stationary multidomain 
pattern in two dimensions the value of e ~ 0.002, so the 
scenarios discussed above are indeed realized. Also, ac- 
cording to the Osipov's criterion j3l[] , the transformation 
of a stationary multidomain pattern into traveling will 
occur when ~ e 2 <C a c . 

For the purpose of completeness let us quote the equa- 
tions obtained in the case of the three-dimensional mul- 
tidomain pattern consisting of spherical domains situated 
on a closed-packing lattice. The derivation follows along 
the same lines as the derivation in two dimensions. In- 
troducing 



3 _ 3C 3 p Q(e sl , Vs ) 

Ph — ' 



(57) 



Airy/la ' 

and rescaling fj and p according to Eq. (p7|), we obtain 



dt 
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(58) 



(59) 



The properties of these equations are essentially the same 
as in two dimensions. 



V. ENTRAINMENT OF THE OSCILLATIONS IN 
A DISORDERED PATTERN 

The multidomain patterns that we studied so far con- 
sisted of domains situated on a regular lattice. We found 
that the motion of the domains occurs on two time scales: 
the short time scale lo^ 1 which is proportional to the pe- 
riod of the oscillations of a single domain and the long 
time scale n. The motion of the domains on the time 
scale Wq 1 is conservative, so the time n is associated 
with the relaxation of the oscillator's energy In two di- 
mensions we found a lot of the situations in which the 
collective domain oscillations should become chaotic, ei- 
ther because of the nonlinear dynamics of the oscillation 
amplitude and phase, or because of the defects. This 
chaotic dynamics, however, is realized on the long time 
scale and is essentially determined by the nonlinear relax- 
ation processes. Here we would like to ask the following 
question. Suppose that instead of the domains of a fixed 
radius on a perfect lattice (hexagonal in two dimensions) , 
we have a random arrangement of the domains with a dis- 
tribution of radii. What will be the dynamics of different 
domains on the short time scale a-vT 1 ? 
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To proceed, we will consider a simplified situation and 
ignore the space dependence of all dynamical quantities 
in the problem. One can think that the domains are con- 
sidered in the system whose size is much smaller than 1. 
Let n(p, t) now be the distribution function of the do- 
main radii p. The distribution n(p, t) must satisfy the 
Liuville equation obtained from Eq. Rescaling ap- 

propriately n, fj, p, and t, and neglecting the terms which 
are not significant on the time scale lOq 1 , from Eq. ( |lS| ) 
and ([22]) we get 



at 



Pq- P d n{p, t)dp, 



dn _dn 

~dt" r] dp~' 



(60) 



(61) 



where pg is a constant which comes from the term 
Q(9 s i,r] s ) in Eq. ( |l8| ) and the distribution n(p,t) is nor- 
malized to unity: 



n(p, t)dp = 1. 



(62) 



The last condition implies the conservation of the num- 
ber of the domains as a function of time. This condition 
is in fact not always satisfied, since we have an absorb- 
ing boundary condition for n at p = and even more 
sophisticated situation at large p because of the possi- 
bility of the domain merging. It is clear, however, that 
this condition will be satisfied if the distribution n has 
finite support at all times and the possibility of merging 
is excluded. 

Let us introduce the quantities 



too = 1, m k (t) = / p k n(p,t)dp, 



(63) 



where k = 1,2, ... ,d. Then from Eqs. (roQ) and 
follows: 



dm^ 
dt 



-kr\m 



fe-i, 



dl l d 
dt =P °- md - 



(64) 



(65) 



Equations ( p4J) and ( |65| ) can be solved exactly. One 
can use Eq. (|64|) with k = to eliminate fj from Eq. 
(p4) with k = 1 and integrate this equation to get TO2 = 
to| + Ci, where C\ is a constant of integration. This 
equation, together with Eq. (|64|) with k — can be 
used to integrate Eq. (64) with k = 2 to get to 3 = 



mf + 3CiTOi + C2, and so forth. Thus, the whole set of 
Eqs. (64) and ( |65| ) can be reduced to a single equation 
for mi. The constants of integration are determined by 
the initial distribution function no(p). If one knows the 
solution of the equation for mi, one then gets 



n(p, t) =n (p- toi (t) + mi (0)) . 



(66) 



_The first conclusion one draws from the solution of Eqs. 
and ([55]) is that the oscillations of different domains 
in a random pattern are always entrained. This is in fact 
obvious right from the start, since the rate of change of 
the radii for each domain is proportional to fj which in 
turn depends integrally on the distribution of radii. How- 
ever, the dynamics of the domains is determined by the 
initial distribution no and may be qualitatively different 
for different initial conditions, even if the parameters of 
the system are the same. 

Let us see what happens when d — 1,2, 3. In the case 
d = 1 the equation for mi is 



d 2 TOi 



+ mi - Pq = 0. 



(67) 



One can see that for any distribution uq (with finite sup- 
port and and assuming that no domain merging occurs) 
the motion is that of a simple harmonic oscillator. The 
energy of the oscillator is determined by the value of fj 
at t — 0. One can see that dynamics of the random pat- 
tern on the time scale uiq 1 is identical to the that of the 
ideally periodic pattern. 

In d = 2 and d = 3 the situation becomes somewhat 
more complicated. In d = 2 we have 



mi 



dt 2 



Ci 



2 
Po 



0. 



(68) 



Here again the dynamics is equivalent to that of an ideal 
hexagonal pattern, but the equilibrium point is shifted. 
The shift is determined by C\ = (p 2 ) — (p) 2 at t = 0. One 
can see from Eq. ( |68"| ) that when G\ > Pq the oscillations 
are not possible. In d — 3 we have 



d 2 m\ 
~d&~ 



3CiTOi + C 2 



Pi 



0, 



(69) 



and the dynamics of the domains may be different from 
that of the domains in an ordered pattern. 

The main conclusion that follows from the analysis 
above is that the disorder, especially small, should not 
significantly affect the dynamics of the pattern and that 
the equations for fj and p should adequately describe the 
behavior of irregular multidomain patterns as well. It is 
also clear that the distortions of the domain shapes will 
not have significant effect either. 



VI. SIMULATIONS OF A CONCRETE MODEL 

In this section we present the result of the numerical 
simulations of the model with a simple cubic nonlinearity: 



Q = e + n-A. 



(70) 



(71) 
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Recently, Muratov and Osipov performed extensive nu- 
merical simulations of this model in two dimensions [fl6| . 
The emphasis of their work was on the instabilities of the 
localized solitary structures (AS) and the evolution of the 
patterns excited by a short localized external stimulus in 
the stable homogeneous system. They were able to con- 
struct a state diagram showing what kinds of patterns 
are realized for different values of a and A for a fixed 
value of e <C 1 from localized initial conditions. Here we 
perform numerical simulations of Eqs. @ and (Q) with 
© and 0) with non-localized initial condition. 

The homogeneous state of the system under consider- 
ation is stable when \A\_> 1/3^3 ~ 0.19. The various 
constants involved are |16|26fl 



2, B = 4, Z 



('1 



c 2 = 0. 



(72) 



Unfortunately, the direct simulations of Eqs. (g) and 
(|J) with e <C 1 in the system of sufficiently large size are 
extremely time-consuming even on a very fast computer. 
In order for the equations for p and fj to be in quantitative 
agreement with the simulations we should have e < 0.01. 
Nevertheless, it is expected to have qualitative agreement 
with the predictions of the preceding sections. All our 
simulations were performed with e = 0.05. 

Our first observation is that if a stable stationary mul- 
tidomain pattern is taken as an initial condition in the 
run with sufficiently small a, the pattern will destabilize 
with respect to the synchronous oscillations of the do- 
mains, and its evolution will depend on how small the 
value of a is and on whether the homogeneous state of 
the system is stable. When a is slightly below the thresh- 
old value and the pattern's period is not very small, the 
oscillations of different domains will synchronize. If the 
value of a is smaller, the domains may start to merge 
during the oscillations, what will result in the formation 
of an irregular pulsating pattern. If a is even smaller 
and the homogeneous state of the system is stable, the 
pattern will collapse into the homogeneous state after 
a few periods of oscillations or transform into a turbu- 
lent pattern. The turbulence here is induced by the self- 
replication of the domains and is qualitatively different 
from the chaotic behaviors discussed above |1(|. If the 
value of a is small enough and the homogeneous state of 
the system is unstable with respect to the uniform self- 
oscillations, the multidomain pattern will collapse into 
the uniform self-oscillations. This is in agreement with 
the predictions of our analysis. 

Besides the synchronous oscillations of the domains, we 
predicted traveling waves of the oscillation phase. These 
traveling waves are indeed realized in the simulations 
(Fig. ||). One has to use the special initial conditions 
to excite a traveling wave. The distributions of 9 should 
consist of the domains whose radii vary smoothly along 
the x-axis, and the distribution of r\ should follow this 
variation with the phase difference of n/2. 

In addition to the plane waves, target waves of the 
oscillation phase are observed. Figure shows a portion 



of a target wave emanating from the top left corner of 
the system. The wave breaks down after a while, as the 
domains begin to merge during the oscillations. As a 
result, an irregular chaotically oscillating pattern forms 
in the system. Notice that eventually, when no merging 
occur any longer, the average size of the domains becomes 
greater than the size of the domains at the beginning. 
This can be interpreted as a self-organized increase of 
the patterns period leading to the stabilization of the 
pattern's oscillations. 

In another simulation in which the value of a was 
smaller, the target pattern broke down not only because 
of merging, but also because of the collapse of smaller 
domains (Fig. |l(J). As a result of the domain collapse, 
patches of the almost homogeneous state formed in the 
system. These patches are then invaded by the domains, 
what leads to the formation of an irregular pulsating pat- 
tern. For smaller values of a this may also lead to the 
formation of autowaves and turbulence induced by self- 
replication of domains Jl6[ . 

In the simulations above we used a hexagonal pattern 
as the initial condition. Let us now see what happens in 
the more realistic situation when the multidomain pat- 
tern is disordered. Figure |Il] shows the evolution of such 
pattern. The initial stationary pattern was in fact gen- 
erated as a result of the Turing instability of the homo- 
geneous state at large a. In the simulation of this figure 
the value of a is small enough, so that the stationary 
multidomain pattern is unstable, and a small disturbance 
was added to 77 in the vicinity of the top left corner of the 
system in order to make the amplitude of the oscillations 
in this region bigger. We see that at early times (the first 
three plates in Fig. |Il|) the oscillations of the domains 
are indeed entrained, despite the differences in the do- 
main sizes and shapes. Notice that in this simulation the 
homogeneous state of the system is unstable with respect 
to the uniform self-oscillations. One can see that here 
the amplitude of the oscillations increases until merging 
and collapse of the domains occur at the top left corner, 
which are then followed by the formation of the uniform 
sclf-oscillations. The uniform self-oscillations eventually 
invade the whole system and destroy the multidomain 
pattern. Yet the uniform self-oscillations and the mul- 
tidomain pattern may coexist for quite a long time. This 
effect was observed in the dynamics of the Turing pattern 
in the CIMA reaction [H| . 

Above we considered the dynamics of the patterns 
consisting of circular hot domains in the cold system. 
Clearly, qualitatively the same dynamics is expected for 
the cold domains in the hot system. In fact, by changing 
9 — > —9, r\ — > —T], and A — ► —A in the concrete sys- 
tem under consideration we will transform hot domains 
into cold and vice versa. A possibility exists for such 
oscillations of the multidomain pattern in which the hot 
domains will transform into cold domains and back. This 



process is illustrated in Fig. 12. There a stationary mul- 
tidomain pattern which formed as a result of the Turing 
instability at A = —0.1 when the hot domains are favor- 
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able was taken as the initial condition. The simulation 
was performed at A = 0.1 when the cold domains are 
favorable. One can see the sequence of transitions from 
the hot domains to cold and back during the oscillations. 
A lot of the domains merge as their size increases, so the 
multidomain pattern quickly transforms into an irregu- 
lar pattern consisting of many disconnected pieces. The 
connectivity of the domains changes with time. After 
t = 1.4 the systems eneters a periodic cycle in which the 
pattern changes periodically from hot to cold. In the run 
with these parameters the pattern eventually collapsed 
into the homogeneous self-oscillations. 

So, in summary, the collapse and merging of the do- 
mains appear to be the major cause for the chaotization 
of the oscillations. The chaos here is due to the underly- 
ing irregularity of the domains themselves. The predicted 
chaotic behavior near the onset of synchronous domain 
oscillations seems to be beyond the capabilities of the 
available computational power. While we were able to 
see steady non-uniform oscillations of a hexagonal mul- 
tidomain pattern of period L p = 2.3 at e = 0.1, a = 0.05, 
and A = —0.1 in the system 40 x 40, the system was yet 
too small to identify these oscillations with chaos. 



VII. CONCLUSION 

Thus, in this paper we studied the dynamics of the 
multidomain patterns of small period in Kf2N systems. 
In order to simplify the problem we took the advantage of 
the smallness of e, the natural small parameter in the con- 
sidered systems |p| 4To| , |l5f . In the limit e — > the dynam- 
ics of the pattern reduces to the free boundary problem 
|jl8| , [26f . This problem may be further simplified by using 
the smallness of the pattern's period. This approach is 
somewhat limited in one dimension, but is always appli- 
cable to higher-dimensional multidomain patterns, if e is 
small enough. In the case of the periodic patterns we 
managed to reduce the number of dynamical variables 
involved to the average domain radius p and the aver- 
age value of r\. The equations for these quantities turned 
out to be universal. The nonlinearities of the original 
reaction-diffusion system enter only via a few numerical 
constants of order 1. In d > 2 the dynamics studied by 
us is in fact the asymptotic limit of the true dynamics as 
e — > 0. Note, however, that the value of e has to be suffi- 
ciently small in order for this approach to be in quantita- 
tive agreement with the actual dynamics. For example, 
for the concrete model studied in Sec. VI this is the case 
only when e < 0.01 (see also Ref. [[l6|). Nevertheless, the 
qualitative agreement is good for 0.01 < e -C 1. 

The analysis performed by us predicts both synchro- 
nization and chaos in the collective domain oscillations of 
the multidomain patterns. There are two types of chaos 
that are realized: the chaos associated with the effects of 
the interplay of the dissipation and dispersion of the non- 
linear waves of the phase of the domain oscillations and 



the chaos associated with the irregularity of the underly- 
ing multidomain pattern. The first kind of chaos is shown 
to correspond to the intermittency and defect turbulence 
in the complex Landau-Ginzburg equation. The second 
kind of chaos is due to the breakdown of the averaged 
dynamics description and is associated with the collapse 
and merging of the domains and the destruction of the 
regular multidomain pattern. Nevertheless, it is shown 
that even in an irregular pattern the domain oscillations 
are synchronized locally, so the chaos is still realized on 
the time and length scales of the dissipation processes. 

Recall that e = yJaDg / D v , where Dg and D v are the 
diffusion coefficients of the activator and the inhibitor, 
respectively. According to its definition, the value of e is 
small as long as a is small and > Dg. However, if 
the values of Dg and D v are of the same order, as is the 
case in typical experiments with the autocatalytic reac- 
tions |37]], this means that a ~ e 2 . At this relationship 
between a and e the stationary patterns are always un- 
stable p|- |l0| , |l^ , |37| , and only autowaves will be realized. 
Our analysis also shows that this will be the case for the 
multidomain patterns. In other words, one should have 
D^ Dg in order for the stationary and more complex 
dynamic patterns (not autowaves) to be feasible in an ex- 
periment. As was emphasized by Epstein and Lengyel, in 
chemical systems, for example, one has to devise special 
methods to make the ratio of the diffusion coefficients 
large p7|| . On the other hand, the value of a is deter- 
mined by the kinetics of the chemical reactions involved 
and can be easily made small. 

It is well known that a solitary pattern (AS) in one 
dimension would destabilize and transform into a pulsat- 
ing pattern when a ~ e, what implies that Dg / D v ~ a 
[|| |l0[|l5 17 3^] . The same condition must be satisfied in 
order for the instability of the homogeneous state of the 
system with respect to the uniform self-oscillations be 
before the Turing instability of the homogeneous state 
|p|- |To| , ^6[ . This is a rather strict requirement for a <C 1. 
On the other hand, as was shown in Sec. IV, the destabi- 
lization of the stationary multidomain patterns in higher 
dimensions occurs when a ~ e 4 / 3 , what implies that in 
this case Dg/D tl ~ a 1 / 2 . This condition requires con- 
siderably smaller difference in the diffusion coefficients of 
the activator and the inhibitor. 

Another interesting implication of our results is that in 
some situations it is possible that the oscillatory behavior 
of the "homogeneous" system is actually the consequence 
of the dynamics of the underlying multidomain pattern. 
Indeed, for the reasonable values of D v ~ 2 x 10 _5 cm 2 s _1 
and ~ 10 s f37| we would have that the size of an in- 
dividual domain is much smaller than L ~ 1.5 x 10~ 2 cm, 
so the domains may actually lie beyond the resolution 
of the experiment. Notice that in the case e <C 1 the 
oscillatory instability of the homogeneous state is always 
accompanied by the Turing instability p[4To| . Numerical 
simulations show that Turing patterns may form as a re- 
sult of the instability of the homogeneous state even if the 
homogeneous state is unstable with respect to the uni- 
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form oscillations and may in fact persist for even smaller 
a Jl^|. So, we suggest that in certain situations it is 
the Turing patterns that exhibit the oscillatory behavior, 
whereas the system's bulk kinetics has relaxation char- 
acter. If this is the case, one would expect to see the 
coexistence of the relaxation and oscillatory kinetics in 
the system, and multiplicity of the oscillation modes in 
a single system with the same parameters. The latter is 
the consequence of the fact that the characteristic time 
and length scales of the oscillations of the multidomain 
patterns strongly depend on the pattern's period which 
is not uniquely determined by the system's parameters. 

The dynamics of the multidomain patterns discussed 
in the present paper was observed by Boissande, Dulos, 
and De Kepper in the experiments on CIMA reaction 
Jl2| |. They were actually able to follow the dynamics 
of the domains and see synchronization, waves of the 
oscillation phase, including spirals on top of the mul- 
tidomain pattern, merging of the domains and the co- 
existence of the domain oscillations and the uniform os- 
cillations of the homogeneous state (the Hopf-holes, as 
they called them). They also emphasized that the waves 
they observed are essentially different from those in the 
Belousov-Zhabotinsky (BZ) reaction. In the CIMA reac- 
tion these are phase waves, as opposed to the autowaves 
in the BZ reaction. This is precisely the conclusion of 
our analysis. 

Synchronously oscillating domains were also observed 
by Rose et al. in the experiments on the catalytic CO 
oxidation on the platinum surface [fblf . The authors sug- 
gested that these oscillations may be explained by the in- 
troduction of global coupling. The results of the present 
paper suggest a more natural explanation for this phe- 
nomenon as dynamics of the multidomain pattern. To 
do this, one has to introduce a diffusion term into the 
equation for the inhibitor in the two-variable reaction- 
diffusion model of this reaction, which is otherwise known 
to have relaxation kinetics |39| . This would also seem to 
explain the transition from a target pattern into a cellular 
structure observed in these experiment. We would also 
like to mention the observation of synchronously pulsat- 
ing cellular flames in the combustion experiments |l3| . 

Haim et al. observed experimentally a single breath- 
ing domain in the FIS reaction in a circular reactor p2[ . 
Their results agree with the conclusions of Sec. IV con- 
cerning the motion of a single domain. Indeed, they see 
the supercritical Hopf bifurcation from the stationary to 
the breathing domain, growing anharmonicity of the os- 
cillations for larger amplitudes of the oscillations, and 
collapse of the domain for yet larger amplitudes. In this 
situation Eq. (50) (with 7r/v3 replaced by 2) might be 
used for the quantitative explanation of these effects. 

The author would like to acknowledge the computa- 
tional support from the Center for Computational Sci- 
ence at Boston University. 
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FIG. 1. The qualitative form of the nullclines of Eqs. (^) 
and (§). 

FIG. 2. The distributions of 6 and r\ in a stationary 
one- dimensional periodic strata. 

FIG. 3. The parameter regions for oscillations and break- 
down of the one-dimensional periodic strata. The upper solid 
line is the instability threshold for the multidomain pattern 
with respect to the oscillations. Below the bottom sold line 
the collective oscillations of the strata break down. 

FIG. 4. The dependence t\ c {j>o) from Eq. ( [52| ) for different 
values of e. 



FIG. 5. The nonlinear friction term f(p) in Eq. (p0[): 
(a) e = 0.00025, n = 0.03; (b) I = 0.0075, n = 0.175; (c) 
e = 0.05, ti = 0.025. Other parameters are Ci = 1, c% = 0. 

FIG. 6. The bifurcation diagram for the two-dimensional 
stationary multidomain pattern for C2 = 0. The numbers 1,2, 
and 3 in the figure correspond to the bifurcations at po = p c i , 
Po = Pc2, and po — p C 3, respectively. The shaded regions are 
those where the instability is not realized (ri e < 0). 

FIG. 7. The sweep of the oscillations for given values of 
Ti (horizontal lines below the solid curve), obtained from the 
solution of Eq. (p6[). The parameters used: e = 0.00025, 
po = 0.25, ci = 1, and C2 = 0. 

FIG. 8. The wave of the phase of the domain oscil- 
lations traveling from left to right. Distributions of the 
activator for different times. The parameters used are: 
e = 0.05, a = 0.02, A = -0.4. The system is 20 x 4. The 
boundary conditions are periodic. 



FIG. 9. Breakdown of a target wave of the oscillation 
phase. Distributions of the activator for different times. The 
parameters used are: e = 0.05, a — 0.019, A = —0.2. The 
system is 20 x 20. The boundary conditions are periodic. 

FIG. 10. Breakdown of a target wave of the oscil- 
lation phase for smaller a. Distributions of the ac- 
tivator for different times. The parameters used are: 
e = 0.05, a = 0.015, A = -0.2. The system is 20 x 20. The 
boundary conditions are neutral. 

FIG. 11. Breakdown of the collective domain oscilla- 
tions into the uniform self-oscillations. Distributions of the 
activator for different times. The parameters used are: 
e = 0.05, a = 0.015, A = -0.1. The system is 20 x 20. The 
boundary conditions are neutral. 

FIG. 12. The oscillations leading to the interconversion 
between the hot and the cold domains. Distributions of 
the activator for different times. The parameters used are: 
e = 0.05, ol = 0.015,4 = 0.1. The system is 20 x 20. The 
boundary conditions are periodic. 
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